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1. Introduction 

Recently, the role of the canonical Boltzmann-Gibbs ensemble as the only useful 
ensemble for quantum statistics has been reevaluated. The need for reevaluation has 
arisen both from experimental and theoretical considerations. Although the canonical 
ensemble works surprisingly well also for systems consisting of only a few particles, 
relativistic ion collisions have provided an arena where the corrections from the finite 
size of the system to the canonical results have become important [Q. A quantitative 
understanding of these collisions is necessary since they are central in the search for 
new phases of strongly interacting matter, such as quark-gluon plasma. In addition, 
methods for a quantitative measurement of the properties of small quantum systems 
could be useful in industrial applications in the future, such as in the design of ever 
smaller electronic components. Similarly, the functioning of certain biological systems 
seems to involve also the quantum behaviour of the system — consider, for instance, the 
recognition of molecules by olfactory receptors p|. 
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The second motivation for non-canonical ensembles comes from the peculiar 
behaviour associated with systems with fractal-like structure. Tsallis statistics 
was developed for the statistical analysis of these cases and we have considered other 
generalized statistics of that kind in a previous work ||^. For "ordinary" systems, 
however, we found out that the so called Gaussian ensemble is the simplest generalization 
of the canonical ensemble which could be used for the approximation of any precanonical 
system — these are systems for which the occupation of high-energy states decreases 
faster than exponentially. For a more complete description of precanonical and Gaussian 
ensembles, see @] and , respectively. Of the non-canonical ensembles, we concentrate 
here only on the Gaussian ensemble with the understanding that the results presented 
can be easily modified to apply for any precanonical ensemble. 

In this work, we develop and prove the existence of a lattice approximation for a 
system of a fixed number of non-relativistic particles which may be fermions, bosons 
or a combination of both. The potential of this system can be very general, the only 
requirement being a fast enough increase at infinity. We do not consider the possibility 
of a creation or annihilation of particles — this would require the use of a quantum 
field theory instead of ordinary quantum mechanics and we would have to give up 
mathematical rigour in the transition. Formulations of the microcanonical ensemble 
using quantum field theory and its perturbative expansion have already been proposed 
by Chaichian and Senda and our approach owes a lot to their work. 

We begin by defining notations and by stating assumptions we need to make on 
the potential in section ^ We then state a theorem which defines and proves the 
convergence of the lattice approximation of traces at a complex temperature and we 
derive a set of lattice observables which can be used for the measurement of energy 
moments. The following section |^ presents a Monte Carlo algorithm for the generation of 
the canonical lattice distribution and we also try to develop an intuitive understanding 
of it. Section ^ contains similar derivations for the Gaussian ensemble, especially it 
contains an essentially canonical algorithm which can be used for the computation of 
the density of states. We conclude with an example where we test the algorithms 
in practice and comment on their strengths and draw-backs. We have also gathered a 
comparison of our methods to the ones presented in the literature along with suggestions 
for improvements in section ^. A good review of the state of the art of using Monte 
Carlo simulations for path- integrals of quantum fluids has been given by Ceperley 0. 

2. Notations and background 

Consider a system of particles living in a (i- dimensional space. In classical physics, this 
system is described by trajectories of these particles and the positions of the particles 
form the relevant degrees of freedom, an n-dimensional space with n = Nd. We shall 
here choose the convention that the positions of particles are denoted by boldface letters 
with a subscript identifying the label of the particle while the classical configuration of 
the system is denoted by regular italics. For example, the n-dimensional configuration 
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X could be written in terms of the N position vectors as x = {xi, . . . ,X]\f). 

The lattice regularization of the quantum effects on this system is defined on a space 
containing L copies of the classical configuration space. The "lattice configurations" are 
thus nL-dimensional vectors and we shall put the number of the copy in parenthesis 
after the configuration symbol. In other words, the position of the particle i in the copy 
k in the lattice configuration x is denoted by Xi{k). The quantum results are recovered 
in the continuum limit, by which we mean here simply taking L to infinity in the lattice 
regularized formulae given below. 

The following discussion applies to massive, non-relativistic particles which are 
bounded by a suitably strong potential and thus can be described by equilibrium 
ensembles. To simplify matters, we also assume that all masses have been scaled away 
by the transformation / ^JrrTi and we use the natural units where the Planck and 

Boltzmann constants are equal to one. 

In summary, we assume that the Hamiltonian of the system is given by 

^ 1 

H = Y^ -Pi"^ + V{xi, ...,Xn) 

4 = 1 

and the precise conditions for the interaction potential are 

(i) V is continuous and bounded from below. 

(ii) / d"i/e-^^(^) < oo for all (3 > 0. 

(iii) V is invariant under permutations of indistinguishable particles. 

With the first condition, the Hamiltonian can be defined in the sense of distributions 
(Theorem X.32 in ^) and the continuity also ensures that the continuum limit will 
be simple — if the original potential is not continuous, it can be replaced by a suitable 
continuous approximation. The second condition requires a sufficiently fast growth 
of the potential at infinity and it gives the precise meaning for what we meant by a 
sufficiently binding potential earlier. The third requirement is more of a consistency 
condition than a restriction. 

The classical partition function of this system is defined by 

^'^^^ - J (27r/5)«/2 ^ 

and clearly ([i|) is equivalent to the assumption that the classical partition function 
converges for all positive temperatures. It is also clear that if V is positive, which could 
be achieved by adding a constant to the potential, then Zc\{l3) decreases monotonically. 

For the sake of simplicity, we shall now assume that there is only one kind of 
indistinguishable particles in the system. The generalization of the results to systems 
with many types of particles should be obvious. 

For each element s in the permutation group Sn we can define a linear 
operator which performs the corresponding particle permutation by s(a:;i, . . . , a^Tv) = 
(a:;s(i), . . . , a3s(jv)). These operators form a unitary representation of the permutation 
group, i.e. all operators are orthogonal, s = 1. The determinant of an operator will 
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also be equal to the parity of the corresponding permutation, denoted here by 
With these notations, the projection operator P± on the physical Hilbert space can be 
defined by 



where (+)* = +1 and the upper sign is used for bosons and the lower for fermions. The 
trace of any trace-class operator A over the physical subspace can also be obtained from 

Trphys(^) = Tr(p±l) = Tr(lP±) . 

Note also that condition (|i|) is equivalent to requiring that V{sx) = V{x) for all 
permutations s. 

3. Definition and convergence of the lattice approximation 

We shall now state a theorem which can be used for getting lattice observables for 
energy moments and which proves the convergence of the continuum limit of the lattice 
approximation. We prove the theorem for the canonical case with a temperature that can 
have an imaginary part — this complex temperature trace is needed for the computation 
of the energy moments in precanonical ensembles as we shall see later. The proof of the 
theorem is given in [Appendix B| . 

Theorem 1 Define T{z) = exp{—z^p^) exp{—zV{x)) on the half-plane Rez > 0. Then 
for all k > and Rez > 0, 

Tr(p±^V^^) = \unJ-l)'^TT(^Pj{z/Lf) . 

All of these functions are analytic on the right half-plane and the convergence is uniform 
on compact subsets of the half-plane. In addition, for any fixed k and < c < c' , these 
functions are uniformly bounded on the closed strip c < Re 2; < c' of the half-plane. 

The proof contains the following explicit representation for the traces on the right 
hand side 

Tr(p±f(;./L)^) = ^ ^(±)^ /■d"^a;ira.a;(L),a;;z), (2) 

where the "lattice kernel function" is defined by 

KL{b,x;z) = exp - -Piib^x) - zVl{x) 

and the "lattice kinetic energy" P^ and the "lattice potential energy" Vl are 

L ^ 

PL{b,x) = -x{k - 1)^, with a;(0) = b, 



2 

k=l 

1 ^ 

vl{x) = -J2nxik)). 

k=l 
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This result is well-known and widely used, it was presented already in the classic work 
by Feynman and Hibbs |p. However, to our knowledge it has not been rigorously proven 
before that the lattice approximation of the complex temperature trace does not spoil 
any part of the result. Note that theorem ^is not a simple consequence of the Trotter 
product formulae |Ty] as these imply only strong convergence for our generators: we 
have to do more work for proving the convergence of the traces. 

The sum in (H) contains A^! terms and thus its computation would be impossible for 
even fairly small values of A^. However, the result of the integral in fact only depends 
on the conjugate class of the permutation and thus the increase of the number of 
terms to compute can be reduced to an exponential e^ increase. To be precise, it is 



shown in [Appendix Q , that 



Ai 



1 ^ 

Trfp±f(z/L)^) = -5^(±1)^-^ ^ 

■ Ai=l A2=0 



xc(A) I d'^^xKLis(x)x{L),x;z), 




N 



(3) 



where S(a) is a representative of the conjugate class defined by the parameters (A/) as 



explained in [Appendix Q and C(a) counts the number of elements in the conjugate class, 
as given by equation [C.2[ . The sum goes over those decreasing sequences of A^ non- 
negative integers whose sum equals A^; in the above equation, 6 is the Kronecker symbol 
which ensures that the last condition holds. 

If all the terms in the sum are relevant, then only about A^ < 30 can be handled by 
the above formula (for A^ = 30 there are 5604 terms in the sum). However, typically only 
a fraction of the terms are relevant and the sum can be used also for a larger number 
of particles. For real temperatures, the estimation of which terms will be relevant can 
be done fairly easily by using lemma ^ given in Appendix B , According to the lemma, 
for all P>0, 



d"^xirL(sx(L),a;;/3) < 



y 



exp 



-^\sy-y\'-PViy) , (4) 



(27r/3)'^/2 

where the right hand side is much easier to evaluate than the left hand side. This also 
gives one more proof of the well-known result that only the identity permutation is 
relevant for dilute gases, for which the typical classical distance between the particles is 
larger than the thermal wavelength a/Stt^. 



4. Lattice operators for energy moments 

We shall next derive a set of lattice operators whose expectation values will converge in 
the continuum limit to the canonical expectation values of powers of the Hamiltonian. 
By the results of section |^, 

Tr(p±^V^^) =2im ^ J](±)^ f d^'^x {-1)'^^Kl{sx{L),x; z), 
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where we have moved the differentiation inside the integral — it is part of the proof in 
Appendix B| that this is allowed. 



By using the Leibniz rule and a bit of algebra we can compute the derivatives, 



d^ 

{-!)'' -^Kl{sx{L),x;z) = Qk{x; L, z)Kl{sx{L), x; z), 



where 



j=0 i=0 ^ ^ ' 

We call Qk the lattice operator for the k:th energy moment. The complicated form 
might be a bit surprising — the lattice action, after all, is simple and it can be guessed 
from the expected continuum limit. We might thus expect to find a lattice Hamiltonian 
equally easily, the guess ''PHl"=I3Vl - ^Pl which goes to "/^f dr[\/(x(r)) - |x2(r)]" 
in the continuum limit, being most natural. 

In fact, the natural guess is not too far away, since by (|^) 

(3Qi{x- L, (3) = (3Vl - ^Pl + ^Ln, (6) 

measures the expectation value of the Hamiltonian normalized by the inverse 
temperature j3. The important difference is that the other energy moments are not 
measured by Q^. Instead, each power needs a separate "renormalization" term, note 
that the extra factors all contain |Ln which diverges in the continuum limit. 

We have given (|^) in a form from where it is apparent that it is the "kinetic energy" 
part of the Hamiltonian, given by Pl, which gets renormalized. This is very natural, 
since the path integration which gives the continuum results goes over continuous but 
non-differentiable paths and thus the derivatives should diverge in the continuum limit. 
In addition to (P) we have 

Q2ix; L, z) = Ql{x- L,z) + \ {\Ln - 2Pl) , (7) 

z 

which can be used for concluding that both the expectation value and variance of the 
lattice operator P^ diverge as \Lnz in the continuum limit. 

5. Monte Carlo algorithm for the generation of the canonical lattice 
distribution 

We saw in section ^ that the term Pl diverges proportionally to the lattice size L in 
the continuum limit. On the other hand, Vl has a finite continuum limit. Consider 
generating the lattice kernel distribution oc exp(— /JV^ — ^Pl) by a straightforward 



Metropolis algorithm |TT|, i.e. make random, local, changes to the lattice configuration 
and keep the changes with the probability exp(— AS*^). Since Pl diverges, for large L 
most changes will be rejected because of the large jump they cause in Pl and the hit-rate 
of the Metropolis algorithm is not very good. But, on the other hand, Pl is a simple 
quadratic function which does not depend on the system, i.e. on the potential or on the 
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nature of the particles. We shall now present an algorithm which takes advantage of 
these simplifications and which behaves better for large L. 

The aim is to generate the canonical lattice distribution for x, 



Ln 
2 



exp 



~-PL{sx{L),x)-f3VL{x) 



where s is any permutation and /3 > 0. Let us first separate the lattice copy containing 
the jump induced by the permutation and define y = x{L). In the remaining n{L — 1) 
integrals we make the following change of variables from x to m, 

k 

x{k) = c(fc) for all 1 < A: < L - 1, (8) 

1=1 

where c{k) = (1 — j^)sy + jy. The Jacobian of this change of variables is one. As earlier, 
for k = and k = L we define also x(0) = sy and x{L) = y. The inverse relation is then 
u{k) = x{k) — x{k — 1) — j^iy — sy) which shows that 



Pl{sx{L),x) = -Y,\<k)-x{k-l) 



k=l 
L-l 



L-1 



1 



\y - sy\ 



k=l k=l 

A choice of other than a straight path c{k) does not lead to as simple a formula, since 
then there would be a cross term linking y and u which is absent for a path with constant 
increments of c. 

For reasons that will become apparent in the next section, we shall use a separate 
notation (3' for the "kinetic temperature" . Thus with the assumption that at this point 
P' = (3, we can express the lattice integral as 



(27r/3')"/^ 



exp 



X 



X 



(27r)(^-iW2 



'-PV{y) 
e 2 L2 exp 



N 



L-l 



i=l k=l 



ew{-^J2l^{x{k))-V{y)]}, 

k=l 



(9a) 

m 

(9c) 



where we have scaled the temperature away from the lattice fluctuations, i.e. made the 
change of variables u = f3' /Lv. This means that in the last line we have used the 
definition 



k 

x{k) = sy + -{y - sy) + 



ft 

1=1 



v{l), for all 1 < < L - 1. 



(10) 



The algorithm consists of generating each of the distributions (pa|-c) separately: 

(i) Compute the classical integral (^) . This is needed for the normalization of the next 
Metropolis step and, by (|), it can be used also for determining if a computation 
of the full term is necessary. 
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(ii) Use e.g. a Metropolis algorithm to generate the classical distribution in ( p^ ) for y. 

(iii) Generate the distribution for v by applying a Metropolis check with normal 
distributed trials for each coordinate of v. This can be done independently of 
the generation of y and, in fact, each classical coordinate can also be generated 
separately. This problem is equivalent to the generation of a self-recurring one- 
dimensional random walk of L steps, the cumulative sums of v corresponding to 
the positions of the walker. Note that the constants in have already been 
chosen so that the result is a probability distribution. 

(iv) Apply a Metropolis check implied by ( p^ for the generated values of y and v. 

In practice, it is best to store the values of x{k) and to use a separate routine 
for the generation of the trials for quantum fluctuations of each particle. The idea 
behind the algorithm is to generate a path of possible quantum fluctuations for a 
classically distributed particle and to test for the acceptance of the whole fluctuation 
path simultaneously. Since Pl does not enter this acceptance test at all, this avoids the 
large rejectance rate from which a local change made in just one copy suffers. But as the 
proposed change is not local, the evaluation of the last Metropolis check is made more 
difficult and time-consuming — this should be more than compensated by the improved 
acceptance rate. 

Let us also briefly comment on the generation of the classical distribution when s is 
not the identity permutation. Then s has at least one cycle of length i > 1 and assume, 
for simplicity, that it is composed of the first particle labels. By the same change of 
variables as before, i.e. by defining yi = J2]=i '^h "we get then 

i=l 1=2 i=2 

This is independent of Wi, i.e. of one of the particle positions. On the other hand, it 
also implies that the distance between any two particles is of the order of \f^. It is 
thus possible to think of each £-cycle as describing a cluster of C. particles which move 
together in the external potential. Similarly, especially for higher temperatures, it is best 
to generate the combinations u from the Gaussian distribution and use the potential 
only for Ui and, of course, for the Metropolis check. 



6. Lattice evaluation of Gaussian traces 



The Gaussian ensemble has been introduced in , where its connections to the canonical 
and the microcanonical ensemble were also illustrated. The ensemble is defined by the 
unnormalized density operator 



Pe{E) 



V2 



: exp 
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It was proven there that for any observable A which satisfies the condition Tr(|y4|e < 
oo for all /? > 0, we have the following integral representation 

TrUp,(E) ) = / ^e^^'^'+^^Trfle-^^) , (11) 
^ ^ i/3-ioo 27ri V / 

where is a positive constant. 

For numerical computations, the proper choice of (3 is very important. Suppose we 

wish to examine the behaviour of the system near the energy scale Eq using the energy 

resolution e < Eq. Kt E = Eq, it was shown in that the best choice for (3 is the 

unique positive solution to the equation 

E^^B^^^H^m^. (12, 

Tr(P±e-'3^) 

In practice, it is usually easier to begin with the temperature 1//5 and then solve Eq 
from equation (0). The natural energy resolution parameter is then e' = j3e and as 
the energy parameter, it will be easiest to use the scaled difference E' = f3{E — Eq). 
Similarly, the moments are most naturally computed as 

P\{H -E- /3£2)fc)|-- = -H)- E')'^)! r- 

We shall derive the Gaussian lattice expressions for this set of moments — all others can, 
of course, be computed from these. 

With the present assumptions, A = P± satisfies the condition for the use of the 
integral representation (0). Therefore, for any real t and /3 > the following holds, 

r./3+ioo 



(13) 



' Jp-ioo 27ri V 
Applying proposition || then shows that both sides can be analytically continued to 
entire functions and that for all k > 0, 

iv(p±(ij - P - Pe')%{E)) = ^ Tr (p±e*(^-^-^^')p,(E) 
r./3+ioo 



i=0 



2m dt'' 



e3^'"'+"-^Tr(P±e-"^ 
t=o V 



' /3— ioo 

The derivative on the right hand side can be computed from the generating function of 
the Hermite polynomials Hk, 



dk 



P2 

dt^ 



It 



3=0 

Parameterizing the integral as z = (3{1 + ia) and using the scaled variables e' and 
E' we have the result 

P^'Tii^hiH ~ E - l3e^)%{E) 

Lfc/2J 



^ ^ 2^j!(fc-2j)! '-'^'^ ' 
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i-oo2vr V / 

Let us next define the function 

g{a;p) = e'"'^^Tr(^P±e~'3(i+^"^^) , 

using which we have for even k, 

h{E,e;P) = {-If/' Re H ^e'"^ {ea)'^e"^^''^' g{a- P), (14) 

and for odd fc, 



h{E, e- P) = im / e'°^(£«) " g{a- P). (15) 

This shows that if we know the function g in some suitably dense set of points, we can 
approximate Ik for any k, e' and E' by multiplying the known values by the weight 
function and then performing a discrete Fourier-transform. 

Thus g{a] P) contains all information even about the arbitrarily high energy end 
of the spectrum and it should not come as a surprise that its computation is not 
easy in general. Fortunately, it has a lattice approximation and there is a relatively 
straightforward Monte Carlo algorithm for the evaluation of the lattice integral. 

By theorem the lattice approximations 

gL{a;P) = e-^^Tr(p±f (/5(1 + i«)/L)^) (16) 

converge to g{a] /3) as L ^ oo and they are all uniformly bounded in a. If we define 
Ik{E, e; P, L) by replacing g by gL in (0) and ([TsD , we can thus rely on the dominated 
convergence theorem and conclude that 

where E' = P{E + Pe' — H) and e' = Pe. Since the Gaussian partition function measures 
directly the density of states, it is of a special interest and, by the above results, it has 
the lattice approximation 

Tiip^UE)) = lim e-^^'''+f'^pio{E',e'-p,L). 
6.1. Monte Carlo algorithm for the Gaussian integrals 

Let us now go through the lattice algorithm for the evaluation oi gi when a and /3 > 
are given. The definition (|T^ translates into the lattice integral 



gL{a-, /5) = ^ E (±)'e'"^'' / d"^^ Kl{sx{L), x; P{1 + ia)) 

r. { 1 „ arctana\l 
X exp yia [PH - pVi + jPl - \Ln j J , (17) 
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where is a short-hand for the combination f3{l + a^) and = Pl{sx{L),x) and 
Vl = Vl{x) as before. In addition, since the kernel of T is obtained by analytic 
continuation, the arcus tangent in the above should be chosen from the primary branch. 



Now we can also finally explain why we wanted to include the term j3H into the 
lattice integration. Suppose, for simplicity, that of the possible classes of permutations 
one dominates over the others — this is natural if the quantum behaviour of the system 
comes from the formation of correlated particle clusters of fixed size for the temperature 
in question (we saw in section ^ that when a permutation is decomposed into cycles, 
each cycle can naturally be identified as such a correlated cluster of particles). Then for 
this permutation and for small a, the last term in parentheses in (|l^ has an almost zero 
expectation value under the distribution in front of it — this follows from a comparison 
with the definition of Qi in @. Since the a-integration is concentrated near zero, 
choosing f]H as the energy origin thus gives the smallest oscillations for the integrand 
in This also shows, that for a fixed lattice size L, it is best to use the value 

evaluated from the same lattice for H, a continuum extrapolation is not necessary or 
desirable. 

The first term in the lattice integral in (0) gives simply a canonical distribution 
which can be generated by the algorithm explained in section ^. The evaluation of the 
remaining exponential is then a simple computation of trigonometric functions. 

Computation of g for all a would then enable the solution of the full microcanonical 
spectrum of the Hamiltonian. As the energy is increased, the dependence of the position 
of the spectral hues depends on ever finer structure of the potential. Since, eventually, a 
computer cannot hold the values of the potential at the accuracy needed, there must be 
a catch somewhere. In our case, this is most easily seen in the prefactor + which 
must the cancelled by the result from the lattice integral — otherwise g would diverge in 
the continuum limit. This means that as a is increased, the oscillations of the lattice 
observable increase rapidly and the maximum value of a that can be computed this way 
is limited by the number of Monte Carlo sweeps feasibly performed and by the resulting 
inaccuracy of the result from the lattice integration. 

7. An example: charged particles in a quadratic potential 

We now want to check how well the algorithm works in practice. We shall consider 
particles living in a three-dimensional world and bound by a harmonic potential. After 
we have checked that the algorithm works for this analytically solvable case, we shall 
add a Coulomb interaction between the particles. In other words, we shall consider the 
potential 



i.e. arctana G (— ^vr, ivr). 
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where gives the number of particles, uj is the binding frequency and q denotes the 
charge of one particle. 

The scaling properties of the potential part of the action are now simple, 

where 

L N ^ L N N ^ 2 

k=l i=l k=l i=l j=l,j^i I ^ -J^ ^1 

By first scaling the temperature (3 away from the kinetic term, and by subsequent 
differentiation we can thus derive a second set of observables for the computation of the 
energy moments. The energy expectation value can now be measured either by using 
Qi in (P) or by using the observable 

= 2Vq + Wc, (18) 

which has been rescaled back to the original lattice variables. Similarly, the second 
moment is given by (|^) or by 

Q2 = Ql + ^i2VQ - ^Vc) 

and the other observables Qk could be computed equally easily. Note that, again, simply 
taking the second power of the observable giving the energy expectation value does not 
yield the correct result. 

The advantage of using both of these two kinds of observables comes from their 
very different dependence on the lattice kinetic energy. For large L the second set of 
observables is better, since it does not contain any cancellation of two large numbers as 
is necessary for the first set of observables. For all values of L, both methods should 
nevertheless yield mutually consistent results and we, in fact, used this for fine-tuning 
and checking of the Monte Carlo algorithm. 

We have presented the results of the Monte Carlo computation of the energy 
moments of a system of two distinguishable three-dimensional particles at two different 
temperatures in Table 0. The results were computed performing 3 ■ 10^ sweeps on 
an AlphaStation 1000 computer which took a few hours of computer time. The error 
estimates for the moments were computed using the bootstrap-method and the number 
of sweeps was enough for obtaining a clean signal for all of the measured six moments. 

In both cases, we used a lattice with L = 8 steps and it is clear that the results are 
very close to the continuum limit. In fact, the main difference to the classical L = 1 
results comes already at the addition of one classical copy for the quantum fluctuations, 
i.e. already at L = 2. As expected, the convergence to the continuum limit is also faster 
for smaller values of (3. 

We have included the (3=1 case as a warning about the possible systematical 
errors induced by the Monte Carlo algorithm. It is clear that the two lattice observables 
at /? = 1 do not yield as mutually consistent results as at /? = 0.1. This is most clearly 
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Table 1. First six moments {I3^{H — {H))'^) of a two-particle system in a three- 
dimensional space. Only the s = id term has been considered, i.e. the particles have 
been assumed to be distinguishable, and the binding frequency is = 1. The results 
were obtained from an eight-step lattice using either Qk as defined in equation (||) or 
Qk as defined in section For the non-charged case, the known continuum value is 
also given. 
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7100(290) 






Q 


6.481(4) 


5.49(3) 


11.4(3) 


123(5) 


790(100) 


8800(2500) 






c.l. 


6.4919 


5.524 


11.95 


127.6 


804.4 


7664 




2 


Q 


8.735(3) 


5.03(2) 


10.8(2) 


110(2) 


647(23) 


6150(430) 






Q 


8.665(4) 


4.69(3) 


9.8(3) 


98(7) 


730(170) 


9400(4000) 


0.1 





Q 


6.000(3) 


5.99(2) 


11.8(2) 


147(2) 


832(27) 


8930(440) 






Q 


6.000(2) 


5.970(12) 


11.87(15) 


143(2) 


886(40) 


9340(790) 






c.l. 


6.0050 


5.995 


12.00 


143.8 


863.4 


8629 




5 


Q 


6.583(3) 


5.82(2) 


11.6(2) 


139(2) 


852(26) 


8690(370) 






Q 


6.585(2) 


5.79(2) 


11.8(2) 


138(3) 


907(65) 


10300(1500) 



seen from the energy expectation value of the charged case when the two results are not 
consistent within their statistical errors when (3 = 1, but they are consistent for /3 = 0.1 
even though we have increased the charge for the second case. 

The source of these systematical errors is in the a priori distribution which is 
used for getting samples for the Metropolis check. If the a priori distribution does not 
concentrate to the correct region of the configuration space, the resulting probability 
distribution in this region will be coarser than suggested by the number of sweeps used 
and the results, likewise, less accurate. 

For instance, the differences in the non-charged case are caused by the use of 
quantum fluctuations as the a priori distribution although we are already in the low- 
temperature region /? — oo where the lowest energy state dominates. This could 
be remedied by using also sweeps with the usual Monte Carlo sampling where the 
binding potential e~^^^ serves as the a priori distribution. The discrepancy in the low- 
temperature charged case is more problematic, since it is caused also by the singularity of 
the Coulomb distribution. The cure in the second case would be replacing the Coulomb 
distribution with some smooth approximation of it, but the effect of this change would 
then have to be analyzed separately. Of course, it is always possible simply to increase 
the number of sweeps until the final probability distribution is smooth enough, but this 
might not always be feasible. 

We also measured the Gaussian partition function, i.e. the density of states, for 
this system. The measurements were performed by computing gL{a;f3) for \a\ < 1.6 
at intervals of 0.05 on a four-step lattice using f3 = 0.1. Samples were then taken from 
the measured values with the assumption of normally distributed errors and a discrete 
Fourier-transformation with the correct weight was performed to each of these samples. 
We have given the stable results of this analysis in figure for two different values for 
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Figure 1. The base ten logarithm of the Gaussian partition function obtained using 
a four-step lattice with lu = 1 and = 0.1 and the two shown values of the energy 
resolution. For the sake of clarity, the results for e = 1 are plotted three orders of 
magnitude smaller than they were measured. The solid and dashed lines represent the 
continuum results for q — 0; the dashed line has been plotted twice, both at its correct 
and reduced location. 



the energy resolution, e' = 1.5 and e' = 0.1. It is interesting that we were able to get 
accurate results also for the second case, for which the effect of the Gaussian weight is 
almost negligible. 

Otherwise, the results behave exactly as expected. The most accurate results are 
obtained near E' = 0, which corresponds to E' ~ 40 in the first case and to ~ 60 
in the second case. The accuracy detoriates as the energy difference increases, faster in 
the case of smaller energy resolution. The addition of the Coulomb term also seems to 
affect the density of states only at the low energy spectrum. 

The non-charged results were also compared to their continuum counterparts which 
were computed using the known partition function 



If we remember that we used only an L = 4 lattice, the results presented in figure |l] 
are surprisingly accurate. When we later checked the difference between the classical 
L = 1 results and the full quantum results at /? = 0.1 for the values of a we used, the 
differences were almost negligible. This seems to imply that, at least for the harmonic 
potential, the classical density of states does give a highly accurate measurement of the 
high energy density of states of the quantum system. 
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In contrast, we were not able to measure the discrete structure of the energy 
spectrum by measuring since this would have needed impractically high values for 
a. We also checked that for detecting the discrete spectrum, the use of large values of 
L is essential. Thus the present method is suitable only for measurement of quantities 
which do not depend directly on the discrete energy spectrum, but only on the density 
of states. 



8. Prologue: improvements and a brief comparison with existing methods 

A number of references came to our attention after the Monte Carlo simulations were 
finished and in this section we shall suggest some improvements motivated by them. 
We first present a better algorithm for the generation of the self-recurrent random walk 



inspired by the work done mostly by Doll and Freeman, see e.g. ||12[. The second 



subsection contains a discussion and references to other possible improvements. 



8.1. Generation of the random walk by Fourier components 

Suppose we make, instead of (H), a change of variables defined by a discrete Fourier 
sine-transform, 

2 ^ 

<k) = - - c(j)] sin(7rf j) , for alll < A; < L - 1. 



The following "orthogonality" relations, which can be found from ||T3[ or computed 
directly by expressing sines as exponentials, are valid for the discrete sine-transform, 

^sin(7r|j) sin(7rf j) = -6kk', 



L 



^ cos(7r|) for k = k' 



cos(7r|)-cos(7r^) 

where we have assumed that both k and k' belong to {1,2, ... ,L — 1}. Using these 
formulae it is easy to prove that the linear transformation which defines the change of 
variables is always invertible with an inverse 

L-1 

x{j) = c{j) + a{k) sm{'K^k) , for all < A; < L, (19) 

k=l 

and that with the choice c{k) = (1 — j^)sy + jy the kinetic term becomes simply 

L-1 



PL = hy-sy\' + J2\aik)\'L'sm\-^). 



k=l 
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Since equation integrates to unity, we can then also deduce the Jacobian of this 
change of variables and conclude that ( |9^ can be replaced by 

L-l , L-l 

(20) 



^n(L-l)^ I r I 1 

^ ^ k=i k=i 



where 



''^"2L2siC2(g)- ^^^^ 

Therefore, we can replace the Monte Carlo step suggested in item (|Iii|) in section 
^ by a generation of the path x{k) via the coefficients a{k) in equation (]TP|). Since 
all a(fc) have a normal distribution, their generation is fast and straightforward. Note 
also that the discrete Fourier transformation need not slow down the algorithm, since 
all necessary trigonometric functions can be computed before entering the Monte Carlo 
loop. 

8.2. Discussion about the Fourier method and further improvements 

Using Fourier-transformation of the path to express quantum statistical path-integrals 
was suggested already by Feynman 0. Later, Doll and Freeman |12[ used a form 



with a cut-off for the number of Fourier modes in a Monte Carlo simulation of 
energies of quantum mechanical systems and they found a performance better than in a 
conventional Metropolis Monte Carlo. However, the convergence of energy expectation 
values was found to be non-monotonic for this "primitive Fourier method" and a 
relatively high number of Fourier modes, fcmax ~ 50, was necessary before the continuum 
value was approached [|14|, |T5 |. 



In our simulations, we did not see any non-monotonicity of the convergence of the 
lattice energy expectation values and, in addition, we found that a relatively coarse 
lattice was sufficient for the computation of nearly continuum results. Emboldened 
by these results we now suggest that using the version of the Fourier transformation 
given in the previous section should solve the problems related to the approach to the 
continuum limit that are manifest in the primitive Fourier method. 

In other words, if the "Matsubara frequencies" used in the primitive Fourier method 
are replaced by those derived in the previous section, i.e. if we replace 
2 2/5' ^ , P' 



^k = TT^ by ai 



{k-nf ' ^ 2L2sin2(ff)' 

then the Fourier-method should converge for smaller values of fcmax = L — \. Similarly, 
the integrals needed in the previous energy estimators should probably be replaced by 
their discretized lattice versions. This idea is not new but it has not been widely 
used, either. 

Since we did not use as complicated a system as in the references, it would now be 
necessary and interesting to repeat the test of the convergence and the comparison of 



computing times for the different methods that were done in [14| and [O], respectively. 
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For this comparison, we would like to point out that the differentiation which we used 
for producing the second set of observables, Qk in section]^, would for a generic potential 
lead to what is essentially the virial estimator defined in |]12|. Note, however, that our 
method also yields similar estimators for all other energy moments, not just for the first 
one. In addition, both of our lattice observables always give the same result on the same 
lattice, whereas the previous ones agree only in the continuum limit. 

We can now also explain why the Fourier method is better than the conventional 
lattice Monte Carlo: using the Fourier coefficients as the a priori distribution gets rid 
of the kinetic term which we suggested was the culprit for the slow convergence of the 
conventional Monte Carlo algorithm for large L. We can also explain why typically no 
estimators are given for energy moments beyond the first one: since each power of the 
kinetic term requires a separate renormalization, guessing the correct lattice observables 
from continuum expressions is difficult and the effect of using wrong constants can 
easily be drastic. In addition, observables for the higher moments cannot be reduced to 
"coordinate-space observables" (i.e. observables living in only one time-slice) as is the 
case for the virial estimator of the energy for distinguishable particles. 

Since the changes we suggest here make the Fourier lattice integrals equivalent to 
the discrete time-step integrals, the methods used for improving the performance of the 
earlier methods should apply equally easily here. One could e.g. employ classical cluster 
algorithms [0 or use smeared or effective potentials |T9|. One interesting application 
which we did not discuss here at all is in the computation of real-time correlation 



functions, for this see e.g. pOl, BT 



Although the algorithm for the computation of Gaussian traces is not practical for 
resolving the discreteness of the energy spectrum, the lattice formulae derived are valid 
for arbitrarily fine energy precision. By following the steps which where used in [^, we 
could derive a lattice formulation with a complicated oscillating kernel which could be 
used for solving also the low-energy discrete spectrum. How far in energy and for how 
complicated systems this kernel can be used in practise needs testing. 



9. Conclusions 



We have presented a lattice Monte Carlo algorithm for the computation of energy 
moments in the canonical Boltzmann-Gibbs ensemble for a system of non-relativistic 
quantum mechanical particles. In fact, none of our results would have changed if 
we were to add a bounded classical observable A{x) to the trace. Thus the present 
method can be used also for the evaluation of microcanonical corrections to the canonical 
expectation values of such observables as is explained in [Q. The overall performance 
of the algorithm, especially for high temperatures, was very good. In addition to the 
algorithm, we presented rigorous limits which can be used for a quantitative estimation 
of when and how the indistinguishability of particles will begin to affect the results. 

We also proposed a lattice formulation which can be used for the evaluation of 
Gaussian expectation values and the Gaussian partition function of these systems. The 
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algorithm we gave for the computation of the lattice integrals was shown to yield correct 
information about the density of states even for quite coarse lattices. If knowledge about 
the discrete energy spectrum is required, then using analytical methods or a modification 
of the algorithm seems to be necessary. The present method has the advantage that 
other non-canonical ensembles in addition to the Gaussian one, as well as many different 
energy resolutions, can be analyzed from the same lattice data. 

In both cases, we saw that the classical ensembles work surprisingly well and we 
noted how the intuitive understanding of the classical case can help in the computation 
of the full quantum results. The main use of the Gaussian lattice methods would likely 
to be in the detection of phase transitions. Although it is necessary to increase the 
number of particles to infinity before a clean signal for a phase transition appears, such 
effects should be visible for a finite number of particles as well. 
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Appendix A. Analytic integrals 

The following result is a straightforward consequence of textbook results, but since it 
is central to all the results presented in the text, we give a brief derivation for it. The 
main content of the theorem is that an analytically parameterized integral is an analytic 
function of the parameter if the integrand is uniformly dominated by an integrable 
function in every compact subset of the parameter space and that in this case it is also 
always possible to take differentiation of the parameter inside the integral. 

Proposition 2 Let Q be an open subset of the complex plane and let X be a measure 
space which is a -finite with respect to a positive measure fi. Consider a function 
F : X X Q ^ <C, which is measurable in the product space and which defines an analytic 
function z F{x,z) for almost every (w.r.t. ^) x. If for every compact subset K of Q 
there is a function qk G L^if^) such that \F{x,z)\ < gxix) for almost every x and for 
every z & K , then the function f : z \^ j d^{x)F{x, z) is well-defined and analytic in 
n. Then also for all z E and positive integers k, -^fiz) = J dfi{x)^F{x, z). 

Proof: f is clearly well-defined for all z E Q. By the dominated convergence theorem 
it is also continuous and then an application of the Morera's theorem shows that it 
is analytic. Using Cauchy's integral formula for the derivatives of /, followed by an 
application of Fubini's theorem then leads to the given expression for the derivatives. 
For more details see [EH]. □ 



When using the proposition here, the cr-finiteness is obvious since we use it only for 
the counting measure on the positive integers and for Lebesque measures on Euclidean 
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spaces. The measurabihty of F in the product space is also trivial since the functions 
we are dealing with will be continuous. 



Appendix B. Proof of the lattice convergence theorem 

We prove here how theorem |I| follows from the assumptions d) and given in section 
The proof consist of four parts: we first show that the k = trace is an analytic 
function which works as a generating function for the other powers. We then derive 
an integral representation for the lattice traces and use this to prove their analyticity. 
Next we show that the lattice traces are uniformly bounded on compact subsets with a 
bound given in terms of the classical partition function. Since the traces converge for 
positive values of the parameter, theorem |I] follows from an application of the Vitali 
convergence theorem and properties of analytic functions. 

However, let us begin by defining the following two functions of a lattice 
configuration x = . . . , x{L)), 

L ^ 

PL^b, x) = — \x{k) — x{k — with a;(0) = 6, and 

k=l 



V,{x) = \Y.^{x{k)). 



k=l 



In the proof and in practical computations the following bound for their exponentials 
will become useful. 



Lemma 3 For any permutation s G Sn and for all L > 1 and (3 > 0, 



(T^x 



(■ 



L 

< 



Ln 
2 



exp 
d"y 



^Pl{sx{L),x) 



(27r/3)"/2 



exp 



-I3Vl{x) 



(B.l) 



Proof: By the relation between arithmetic and geometric means, we always have 
exp[— ^^^-^ /?V^(x(A;))] < 7; X]fc=i 6^P[~/^^(^(^))]- Applying this and the rules of 
Gaussian integrals then yields 



d' 



X 



Ln 
\ 2 



) 



exp 



1 

k=l 



-Pl{sx{L),x) 



I3Vl{x) 



\2TrP 



exp 



(TycTxe 



L 



2'k(3{L - k) 



exp 



\2Tx(3k 
L 



exp 



2(3k' 



\sx 



\y 



x\ 



2(3{L - k) 

But s is an orthogonal transformation and thus \sx — = \x — s^y\'^. Substituting 
this into the previous formula and computing the resulting Gaussian integral over x will 
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then show that, in fact, all the terms in the sum are equal and the result is precisely the 
middle term in equation ( [B.l| ). This proves the first inequality. The second inequality 
is obvious from the definition of Zd, equation ([l|). □ 

Let us begin the proof of the theorem by defining for all Kez > 
g{z) = Tr(P±e-^^). 

It has been proven in |^ that for the potentials we are interested in the operator e~^^ 
is trace-class on the right half-plane. Since P± is bounded, the trace converges and can 
be given in terms of the eigenvectors ipj and the corresponding eigenvalues Ej of the 
Hamiltonian as the (numerable) sum 



Suppose next that Re^; > c > 0. Since then \{iljj\P±iljj)e ''^^\ < e and X^j-e = 
rjj^g-cH ^ ^g^^ apply proposition || and conclude that g is analytic on the right 

half-plane and that 

9^'\^) = = (-l)'=Tr(p±ifV^^) . (B.2) 

j 

The operator defining the lattice trace, T{z), is for all Re 2; > a Hilbert-Schmidt 
operator with a kernel 

Therefore, P±T{z) is also Hilbert-Schmidt and it has a kernel 
^5:(±)-/f'-)(.a,6). 

Thus for any L>2 the function fiiz) = Ti (^P±T{z/L)^^ is well-defined and it has an 
integral representation 

1 r T -1 

/^(^) = m5^^^^ /^"^K2^) ' ^^p[--Pl{sx{L),x)-zV,{x)]. (B.3) 
■ sgSn J ■nz z 

We shall use equation (p.3|) to define also fi{z). 

Suppose next that Re^; > c > 0. Then the integrand in (|B.3|) is bounded by 

Ln _ 

(L/27rc)~ exp(— cVi(a;)) which is integrable by assumption. Thus the conditions for 
proposition |^ are satisfied and each fiiz) is an analytic function on the right half plane 
and all of its derivatives can be computed by a differentiation inside the integral. 

We shall now first finish proving the theorem for the /c = case, from which the 
A; > case follows quite easily. The uniform boundedness on compact subsets will follow 
from the bound 

, „ I 1 for L even 

I/lWI <e-^^^^^Zci(aLRe2), whereaL = <^ . . (B.4) 



for L odd 
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and bi = (Lmod2)Vmin/-^ with Vmin = i^ixVi^). Note that for even L we have the 
simple bound |/l(^)| < ZdiRez). 

By Lemma IV. 1 of for any bounded operator A and Hilbert-Schmidt operator 
f , I Tr(lf2^)| < ||1|| Tr(ftf)^. Since ||P±|| = 1 and \\f{z)\\ < e"^''^^--, we thus have 

|/l(^)| < e-^^^'=^Tr(f (aiRez/£)0, where i = [L/2\ . 

Choosing s = id in the left hand side of equation ( |B.3|) yields an integral representation 
for the trace in the above equation. By lemma then Tr(T(/5/£)^) < Zci{P) and we 
have finished proving the bound |B.4| . 



We next need the result that for all /5 > 0, 

lim /i(/5)=Tr(p±e-'5^). 

The proof of this result can be done essentially identically to the proof of Theorem III. 4 
in which contains the above result for the case P = 1. We do not reproduce the 
proof here. 

We have now shown that the sequence fi consists of functions analytic on the right 
half-plane, which converge on the real axis and which are uniformly bounded on every 
compact subset of the half-plane. By the Vitali convergence theorem, the sequence then 
converges on the whole half-plane, the convergence is uniform on compact subsets and 
the limit function is analytic. Since the limit function is equal to g on the real axis, and 
g is analytic, it follows that the limit function is given by g on the whole half-plane. 
Also, since the convergence is uniform on compact subsets, then all derivatives converge 
as expected, f^\z) g^^\z). By equation ( [B.2|) we have now completed the proof of 
the first half of theorem |l|. 

Suppose then that < c < c' and z satisfies c < Re z < c'. By the Cauchy estimates 
for derivatives, then for all < t < 1, 

C (l-t)c<Kew<c'+tc 

and it is thus enough to prove the uniform boundedness of fiiz). However, since Zd 
is now a continuous function, this is an obvious consequence of the inequality ( [B.4|) . 
Especially, if V is positive then Z^i is decreasing and we have for all even L, < t < 1 
and Rez >c, 

l/fWI<^^cl((l-t)c). 

Appendix C. Simplification of the permutation sum 

Let r be a permutation and make a change of variables from x to ?/ to the lattice integral 
(0) as defined by x{k) = ry{k) for all A; = 1, . . . , L. The Jacobian of this transformation 
is one and, by assumption, Vl{x) = Vl(?/). The transformation r is orthogonal, and 
thus \x{k) — x{k — 1)P = \y{k) — y{k — l)p for all k with the exception oi k = 1 for 
which we have \x{l) — sx{L)\'^ = \y{l)—r"^sry{L)\'^. Therefore, the lattice kernel satisfies 
KLisx{L), x; z) = KLiir-^sr)y{L),y; z). 
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Since r was arbitrary, we can now conclude that the result of the lattice integral 
depends only on the conjugate class of the permutation. We shall thus need a way 
of parameterizing the conjugate classes and of choosing an easy representative element 
from each class. The solution is well-known and we present the results as they are given 
in 0. 



Every permutation can be decomposed into independent cyclic permutations, 
cycles, and its conjugate class is determined by the number ui of cycles of length /. 
Conversely, a collection of numbers ui, I = 1, . . . , N, defines a conjugate class provided 
it satisfies the consistency condition Ylf=i ~ ^ ■ -^^^ each class, we can then choose a 
representative element as the permutation which performs the cyclic permutations for 
consecutive elements and puts the longest cycles first, e.g. if = 7 and z/3 = 1, z/2 = 1, 
Pi = 2, the representative permutation would be (1,2,3,4,5,6,7) — > (3,1,2,5,4,6,7). 
The number of distinct permutations in each class can be computed from ui by the 
formula 

^' = rr [ r ^ - ^T=\ J '^AiM.] 

Each cycle with an even number of elements has an odd parity and a cycle with an odd 
number of elements has an even parity. Thus the parity of the permutations in the class 
is given by the formula 

For our purposes, there exists a better alternative parameterization. From vi define 
the new parameters A;, / = 1, . . . , A^, via 

N 

j=i 

The inverse relation is clearly given by ui = Xi — A/+i for / = 1, . . . , A^, where it is 
understood that A/ = whenever / > A^. The consistency condition for A; is now simply 

N 

^ A, = AT and Ai > A2 > • ■ ■ > Aiv > 

1=1 

and thus the generation of the conjugate classes reduces to the generation of the 
partitions of A^ into A^ ordered non-negative integers. After this is done, we can compute 
z// and choose the representative permutation, denoted from now on by S(a), as was 
explained in the previous paragraphs. In addition, the number of permutations in a 
class is given by 

AM 

C(A) = (C.2) 

and, by ( |C.1| ), the parity of the class is given by (— l)^^'^^ 
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